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1 INTRODUCTION 



1.1 Motivation 



The Modified Newtonian Dynamics (MOND) -Milgrom 
(1983a, b,c)- is a successful alternative framework for in- 
terpreting the galactic rotation curves and the empirical 
Tully-Fisher relation without relying on dark matter halos 
(see Sanders & McGaugh 2002 for a review). At the non- 
relativistic level, the best modified-gravity formulation of 
MOND is the modified Poisson equation originally proposed 
by Bekenstein & Milgrom (1984), 



fflo 



where p is the density of ordinary (baryonic) matter, U is 
the gravitational potential, g = VC7 is the gravitational 



* E-mail: blanchet@iap.fr 

f E-mail: jerome.novak@obspm.fr 



ABSTRACT 

The Modified Newtonian Dynamics (MOND) has been formulated as a modification 
of the Poisson equation for the Newtonian gravitational field. This theory generically 
predicts a violation of the strong version of the equivalence principle, and as a result 
the gravitational dynamics of a system depends on the external gravitational field in 
which the system is embedded. This so-called external field effect has been recently 
shown to imply the existence of an anomalous quadrupolar correction, along the di- 
rection of the external galactic field, in the gravitational potential felt by planets in 
the Solar System. In this paper we confirm the existence of this effect by a numerical 
integration of the MOND equation in the presence of an external field, and compute 
the secular precession of the perihelion of planets induced by this effect. We find 
that the precession effect is rather large for outer gaseous planets, and in the case of 
Saturn is comparable to published residuals of precession obtained by Saturn range 
tracking data. The effect is much smaller for inner planets, but in the case of the 
Earth it appears to be in conflict for most of the MOND functions /i(y) with the very 
good constraint on the perihelion precession obtained from Jupiter VLBI data. The 
MOND functions that are compatible with this constraint appear to have a very rapid 
transition from the MONDian regime to the Newtonian one. 

Key words: planets and satellites:general, galaxies:kinematics and dynamics, meth- 
ods:numerical, dark matter 



field and g — \g\ its ordinary Euclidean norm. The modi- 
fication of the Poisson equation is encoded in the MOND 
function n(y) of the single argument y = g/ao, where 
ao = 1.2 x 10~ 10 m/s 2 denotes the MOND constant acceler- 
ation scale. The MOND function interpolates between the 
MOND regime corresponding to weak gravitational fields 
y = g/ao <C 1, for which it behaves as fj,(y) = y + o(y), and 
the Newtonian strong-field regime y ^> 1, where /i reduces 
to 1 so that we recover the usual Newtonian gravity. 

Relativistic extensions of MOND modifying general rel- 
ativity have been proposed (Bekenstein 2004; Sanders 2005) 
and extensively studied (see Bekenstein 2005 for a review). 
Alternatively the MOND equation (1) can be reinterpreted 
as a modification of dark matter rather than gravity by 
invoking a mechanism of "gravitational polarisation" — a 
gravitational analogue of the electric polarisation (Blanchet 
2007) ; this yields to the concept of dipolar dark matter which 
has been formulated as a relativistic model in standard gen- 
eral relativity by Blanchet & Le Tiec (2008, 2009). 

An important consequence of the non-linearity of 



= -47rGp , 



(1) 
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Eq. (1) in the MOND regime, is that the gravitational dy- 
namics of a system is influenced (besides the well-known 
tidal force) by the external gravitational environment in 
which the system is embedded. This is known as the exter- 
nal field effect (EFE), which has non-trivial implications for 
non-isolated gravitating systems. The EFE was conjectured 
to explain the dynamics of open star clusters in our galaxy 
(Milgrom 1983a, b,c), since they do not show evidence of 
dark matter despite the involved weak internal accelerations 
(i.e. below ao). The EFE effect shows that the dynamics of 
these systems should actually be Newtonian as a result of 
their immersion in the gravitational field of the Milky Way. 
The EFE is a rigorous prediction of the Bekenstein-Milgrom 
equation (1), and is best exemplified by the asymptotic be- 
haviour of the solution of (1) far from a localised matter 
distribution (say, the Solar System) , in the presence of a con- 
stant external gravitational field g c (the field of the Milky 
Way). At large distances r = |x| — > oo we have (Bekenstein 
& Milgrom 1984) 



by 



U = g c ■ x + 



GM/fi c 



ry/l + A c sin 2 9 



+ o 



(2) 



where M is the mass of the localised matter distribution, 
where 9 is the azimuthal angle from the direction of the 
external field g c (see also Fig. 1), and where we denote 
u c = n(y c ) and A c = y c n' c /n c , with y c = g c /a and fi' c = 
d/i(?/ c )/d?/ . In the presence of the external field, the MOND 
internal potential u = U — g e ■ x shows a Newtonian- like fall- 
off ~ r _1 at large distances but with an effective gravita- 
tional constant G/ fie- 1 However, contrary to the Newtonian 
case, it exhibits a non-spherical deformation along the direc- 
tion of the external field. The fact that the external field g c 
does not disappear from the internal dynamics can be inter- 
preted as a violation of the strong version of the equivalence 
principle. 2 For the reader's convenience we derive the result 
(2) in the Appendix A. 

In a recent paper, Milgrom (2009) has shown that the 
imprint of the external galactic field g e on the Solar System 
(due to a violation of the strong equivalence principle) shows 
up not only asymptotically, but also in the inner regions of 
the system, where it may have implications for the motion 
of planets. This is somewhat unexpected because gravity is 
strong there (we have g ao) and the dynamics should 
be Newtonian. However, because of the special properties of 
the equation (1), the solution will be given by some modified 
Poisson-type integral, and the dynamics in the strong-field 
region will be affected by the anomalous behaviour in the 
asymptotic weak-field region. Thus the results apply only to 
the nonlinear Poisson equation (1), not to modified inertia 
formulations of MOND. The anomaly expresses itself as an 
anomalous quadrupolar contribution to the MONDian inter- 
nal field u, as compared to the Newtonian potential, given 



1 Recall that in the absence of the external field the MOND po- 
tential behaves like U ~ — \/GMciq lnr, showing that there is 
no escape velocity from an isolated system (Famaey et al. 2007). 
However since no object is truly isolated the asymptotic behaviour 
of the potential is always given by (2), in the approximation where 
the external field is constant. 

2 However the weak version of the equivalence principle, that all 
test particles in the MOND gravitational field have universal ac- 
celeration a = g, remains satisfied (Bekenstein & Milgrom 1984). 



8u = -x l x j Qij , 



(3) 



where x % is the distance from the Solar System barycentre, 
and Qij is a trace-free quadrupole moment of the type 



Q^j = Q 2 (eV - i^j) , 



(4) 



with Q2 being the quadrupole coefficient and e = (e 1 ) the 
preferred direction of the galactic centre (i.e. e = g e /g e ). 
The anomalous term (3) is a harmonic solution of the 
Laplace equation, i.e. ASu — 0. 

The radial dependence of this anomaly is oc r 2 and can 
thus be separated from a quadrupolar deformation due to 
the Sun's oblateness which decreases like oc r™ 3 . This result 
is valid whenever r is much less than the MOND transition 
distance for the Solar System, defined by r = ^/ GM /ao 
with M the mass of the Sun and ao the MOND accelera- 
tion scale. This radius corresponds to the transition region 
where the Newtonian acceleration becomes of the order of 
the MOND acceleration ao and therefore, MOND effects be- 
come dominant. We have ro — 7100 AU so the effect (3)-(4) 
is valid in a large volume around the Sun including all the 
planets (recall that Neptune's orbit is at 30 AU). 

In addition to (3) we have also the usual MOND ef- 
fect. In the special case of spherical symmetry, the equation 
(1) reduces to fi(g/ao)g = <?n, with the ordinary Newto- 
nian gravitational field. For a MOND function behaving like 



Mi/) 



«2/ 



when y — > 00 (i.e. in the region r <IC r ), 



the Newtonian field is modified by 5gw ~ eao(r/ro) 27 ~ 2 . For 
7=1 this corresponds to a constant acceleration similar to 
the Pioneer anomaly, and for 7^2 the effect is very small 
and the motion of planets is almost unaffected. This effect 
is independent of the presence of the external field, and is 
spherical, so can be distinguished from the quadrupolar de- 
formation (3) induced by the external field; we neglect this 
effect here. 



1.2 Summary 

In the present paper we shall confirm the existence of the ef- 
fect (3)-(4) and shall derive it by an independent numerical 
integration of the equation (1) using an elliptic solver orig- 
inally built for numerical relativity purposes. 3 The magni- 
tude of the quadrupole coefficient Q2 depends on the dimen- 
sionless ratio between the external field g c and the MOND 
acceleration a , say 



5c 
ao 



(5) 



and on the particular MOND interpolating function « in 
use. For the Milky Way field at the level of the Sun we have 
g c ~ 1.9 x 10~ 10 m/s 2 which happens to be slightly above 
the MOND scale, i.e. 77 ~ 1.6. Our calculation implemented 
for various standard choices for the /^-function gives 



2.1 x HT 27 s" 2 < Q 2 < 4.1 x HT 26 s~ 2 . 



(6) 



3 The code is based on the library LORENE, available from the 
website http : //www . lorene . obspm . f r . 
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The quadmpolc coefficient Q2 is found to be positive corre- 
sponding to a prolate elongation along the quadrupolar axis. 
We shall verify numerically that Q2 is indeed constant over 
a large portion of the inner Solar System, and starts decreas- 
ing when the distance becomes comparable to the MOND 
transition scale ro, at a few thousands AU from the Sun. 
We shall also check that Q2 scales approximately like the 
inverse square root of the central mass, in agreement with 
dimensional analysis which shows that we should approxi- 
mately have Q2 = 92*10/7*0, where 92(77) is a dimensionless 
quadrupole coefficient depending on the ratio rj between the 
external field and the MOND acceleration. We find that for 
rj ~ 1.6, q2 ranges over 0.02 < q2 < 0.36, in good agreement 
with the values found by Milgrom (2009). 

We then study some consequences of the anomalous 
quadrupole moment for the motion of planets, notably the 
secular advances of planetary perihelia. Applying standard 
perturbation equations of celestial mechanics we obtain the 
anomalous precession rate 

To simplify, we have assumed that the direction of the galac- 
tic centre lies in the plane of the ecliptic (this is only approx- 
imately correct) and that all planets move in this plane. We 
defined the origin of the precession angle Co to be the di- 
rection of the galactic centre. The orbital frequency of the 
planet is n — 2n/P (P is the orbital period), and a, e and 
u = u + fl denote the standard orbital elements. The effect 
seems to be the most interesting for Saturn for which we 
find 

0.3 mas/cy < A 2 < 5.8 mas/cy. (8) 

Those values are within the range of published residuals 
for the precession of Saturn as obtained from global fits of 
the Solar System dynamics in Pitjcva (2005) and Fienga 
et al. (2009). However we shall find that in the case of the 
Earth the predicted perihelion precession may be greatly 
constrained by the best estimates of postfit residuals ob- 
tained thanks to the Jupiter VLBI data using the INPOP 
planetary ephemerides (Fienga et al. 2009). This is turn re- 
duces the possible choices for the MOND function n(y) to 
those that exhibit a sharp transition from the Newtonian to 
MONDian regime. 

To conclude, we present very accurate numerical com- 
putations performed within the well-defined theory (1), and 
compare the results with recent observations of the motion of 
planets in the Solar System. We find that the observational 
constraints are rather strong, and may even conflict with 
some of the predictions. Thus our results provide motivation 
for investigating more systematically possible anomalous ef- 
fects in the Solar System predicted by alternative theories. 
In particular, the external field effect associated with the 
preferred direction of the galactic centre could be seen as 
a typical one arising in generic attempts to modify gravity 
with motivation coming from dark matter and/or dark en- 
ergy. Indeed we expect that generic departures from general 
relativity will violate the strong equivalence principle (Will 
1993). 

On the other hand let us cautiously remark that MOND 
and more sophisticated theories such as TeVeS (Bekenstein 
2004) , which are intended to describe the weak field regime 



of gravity (below ao), may not be extrapolated without mod- 
ification to the strong field of the Solar System. For in- 
stance it has been argued by Famaey & Binney (2005) that 
a MOND interpolating function \i which performs well at 
fitting the rotation curves of galaxies is given by \i\ defined 
by (34a) below. However this function has a rather slow 
transition to the Newtonian regime, given by fii ~ 1 — 
when y = g/ao — ¥ 00, which is already excluded by Solar 
System observations (Sereno & Jetzer 2006). Indeed such 
slow fall-off —y 1 predicts a constant supplementary accel- 
eration directed toward the Sun 5<7n = ao (i.e. a "Pioneer" 
effect), which is ruled out because not seen from the mo- 
tion of planets. Thus it could be that the transition between 
MOND and the Newtonian regime is more complicated than 
is modelled in Eq. (1). This is also true for the modified dark 
matter model (Blanchet & Le Tiec 2008, 2009) which may 
only give an effective description valid in the weak field limit 
and cannot be extrapolated as it stands to the Solar System. 
While looking at MOND-like effects in the Solar System we 
should keep the previous proviso in mind. The potential con- 
flict we find here with the Solar System dynamics (notably 
with the Earth orbital precession) may not necessarily inval- 
idate those theories if they are not "fundamental" theories 
but rather "phenomenological" models. 

The plan of this paper is the following. In Sec. 2, we 
derive an expression of the solution of the MOND equation 
near the origin in terms of multipole moments. Section 3 is 
devoted to numerical techniques and results, with a descrip- 
tion of the particular formulation and assumptions used for 
the numerical integration in Sec. 3.1, of the various tests 
passed by the code (Sec. 3.2), and the numerical results and 
values for the first multipoles of the potential in Sec. 3.3. 
Section 4 details the consequences of this modified gravita- 
tional potential on the orbits of Solar-system planets, with 
first the derivation of the perturbation equations in Sec. 4.1 
and numerical values for the planetary precessions (Sec. 4.2). 
Finally, Section 5 summarises the results and gives some 
concluding remarks. 



2 MULTIPOLAR EXPANSION OF THE 
MONDIAN POTENTIAL 

In this paper we shall solve the modified Poisson equation (1) 
with the boundary condition given by the constant external 
gravitational field g E (defining a preferred spatial direction 
denoted e = g e /g e ), i.e. that the gravitational field g — V(7 
asymptotes to lim, — >00 g = g e . The external field g e should 
consistently obey a MOND equation, but here we shall sim- 
ply need to assume that g e is constant over the entire Solar 
System. The MONDian physicist measures from the motion 
of planets relatively to the Sun the internal gravitational 
potential u defined by 

u = U - g a ■ x , (9) 

which is such that lim^oo u — 0. Contrary to what happens 
in the Newtonian case, the external field g a does not disap- 
pear from the gravitational field equation (1) and we want to 
investigate numerically its effect. The anomaly detected by a 
Newtonian physicist with respect to the MONDian physicist 
is the difference of internal potentials, 

5u = 11 — , (10) 
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where «n denotes the ordinary Newtonian potential gener- 
ated by the same ordinary matter distribution p, and thus 
solution of the Poisson equation Amn = — 4-7rGp with the 
boundary condition that lim,.-,.^ «n = 0. Hence un is given 
by the standard Poisson integral 

u N (x,«) = G j T^^7| P( x '>*)- (11) 

A short calculation shows that the anomaly obeys the or- 
dinary Poisson equation A8u — — 47rGp p dm, where p p dm is 
the density of "phantom dark matter" defined by 

p pdm = • ( X Vf/) , (12) 

where we denote \ = A* — 1- The phantom dark matter rep- 
resents the mass density that a Newtonian physicist would 
attribute to dark matter. In the model by Blanchet & Le 
Tiec (2008, 2009) the phantom dark matter is interpreted 
as the density of polarisation of some dipolar dark matter 
medium and the coefficient \ represents the "gravitational 
susceptibility" of this dark matter medium. 

The Poisson equation A8u — — 47rGp p dm is to be solved 
with the boundary condition that linv^oo 8u = 0; hence the 
solution is given by the Poisson integral 

f d 3 x' 

<5u(x, t) = G J | x _ x ,| Ppdm(x', t) . (13) 

We emphasise here that, contrary to the Newtonian (linear) 
case, the knowledge of the matter density distribution does 
not allow to obtain any analytic solution for the potential; 
however, we can still infer the structure of the multipolar 
expansion near the origin, and the moments will be com- 
puted numerically. We can check that the phantom dark 
matter behaves like r -3 when r — > oo, so the integral (13) 
is perfectly convergent. 

We want to discuss the influence of the external galactic 
field in the inner part of the Solar System where the gravita- 
tional field is strong (g 3> ao); thus p tends to one there, so 
X tends to zero. For the discussion in this section we adopt 
the extreme case where x is exactly zero in a neighbourhood 
of the origin, say for r ^ e, so that there is no phantom 
dark matter for r < e; in later sections devoted to the full 
numerical integration we shall still make this assumption 
by posing \ = inside the Sun (in particular we shall al- 
ways neglect the small MOND effect at the centre of the 
Sun where gravity is vanishingly small). If p p dm = when 
r ^ e we can directly obtain the multipolar expansion of the 
anomalous term (13) about the origin by Taylor expanding 
the integrand when r = |x| — > 0. In this way we obtain 

+°° f-V 

S u = ^2^J-x l Ql, (14) 
where the multipole moments near the origin are given by 4 

Ql = g^ d 3 xp pdm a L (jj , (15) 

4 Our notation is as follows: L = i\ - ■ - ii denotes a multi-index 
composed of I multipolar spatial indices i\ , • • ■ , i\ (ranging from 1 
to 3); 8l = di 1 ■ ■ ■ di, is the product of I partial derivatives di = 
d/dx i ; x L = x 11 ■ ■ ■ x' 1 is the product of I spatial positions x i ; 
similarly n L = n n • • • n l > = x L /r is the product of I unit vectors 
n 1 = x l /r; the symmetric-trace-free (STF) projection is indicated 
with a hat, for instance x L = STF[x L ], and similarly for n L and 



Because the integration in (15) is limited to the domain 
r > e and 9z,(l/r) is symmetric-trace-free (STF) there [in- 
deed A(l/r) = 0], we deduce that the multipole moments 
Ql themselves are STF. This can also be immediately in- 
ferred from the fact that ASu = when r ^ e, hence the 
multipole expansion (14) must be a homogeneous solution 
of the Laplace equation which is regular at the origin, and 
is therefore necessarily made solely of STF tensors of type 
x L . Hence we can replace x L in (14) by its STF projection 
x L . It is clear from the formula (15) that the MONDian 
gravitational field (for r > ro) can influence the near-zone 
expansion of the field when r — > 0. 

With the expression (12) for the phantom dark matter 
and the MOND equation (1), we can further transform Ql 
as 

Q L = -±f^d 3 x[4nGp + Au]d L fy . (16) 

Approximating the central matter distribution as being 
spherically symmetric (i.e. ignoring the "back-reaction" of 
the non-spherical anomalous field 5u on the matter which 
generates the field), we see that the first term is non-zero 
only in the monopolar case I = where it reduces in the 
limit £ ~ > to minus the Newtonian potential evaluated at 
the origin. On the other hand the second term in (16) can 
be transformed as a surface integral. We find 

OL = -u N (0)5(,„+i-^ r °° <fSi [ud iL (^j-diUd L 

(17) 

Our notation means that the surface integral is composed of 
two contributions, an inner one at r = e (denoted Q L ) and 
an outer one at infinity r = oo (say Q2 3 ). In Eq. (17) we 
implicitly assume that the limit e — > is to be taken after 
evaluating the inner surface integral. 

Inserting U — u + g c ■ x into the surface integral at 
infinity we find that it contributes only to the dipolar 
term and reduces to the external galactic field; thus Q°° — gl 
with all other Qff ' s being zero. On the other hand the inner 
surface integral Q L is found to have a well-defined limit when 
e — > given by Q L = (— ) l (8lU)(0), where we recall that 
8l denotes the STF part of 8l = di 1 • • ■ di t . We then find 
that the galactic field in U = u + g c ■ x cancels the dipole 
term in the surface integral at infinity, so that the result is 

Ql = -M N (0)<5i, + (-) l (dLu){0) . (18) 

The internal MONDian potential u admits therefore the fol- 
lowing STF multipolar expansion (equivalent to a near-zone 
expansion when r — > 0) 

u = u^-u^{0) + ^-x L (dLu){Q). (19) 

In Appendix B we derive an alternative proof of this result. 
At this stage we have elucidated the structure of the multi- 
pole expansion of the anomaly 8u near the origin, but still 
we need to resort to a numerical integration of the non-linear 
MOND equation in order to obtain quantitative values for 

&L- The decomposition of 8l in terms of STF components 8l is 
given by (B2)-(B3) below. In the case of summed-up (dummy) 
multi-indices L, we do not write the I summations from 1 to 3 
over their indices. 
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Figure 1. Setting of the spherically symmetric star and the 
asymptotic galactic field g c , using spherical coordinates {r, 8, ip}. 



the multipole moments. Section 3 will be devoted to this 
task. 

Finally we show that the dipole moment Qi (with I — 1) 
is actually zero (Milgrom 2009) . This follows from the weak 
equivalence principle satisfied by the Bekenstein-Milgrom 
theory. As a consequence, the total acceleration of the cen- 
tre of mass of the Solar System in the galactic external field 
g e does not depend on the Solar System's internal dynamics 
and is simply given by g e (for a detailed proof see Bekenstein 
& Milgrom 1984). The centre of mass of the Solar System 
does not differ much from that of the Sun, so we deduce 
that the "self-acceleration" of the Sun, i.e. the total accel- 
eration due to the internal potential u when integrated over 
the whole volume of the Sun, must be (approximately) zero: 

j d 3 xp3iU = 0. (20) 

We have obviously the same result for the Newtonian poten- 
tial mn in Newtonian gravity so the same will be true for the 
anomalous field Su — u — wn- In other words the phantom 
dark matter which is the source for the anomaly exerts no 
net force on the Sun and we get, for a spherically symmetric 
star, 

Qr = ~]g y d 3 xpcWu = . (21) 

In Sec. 3.3 we shall numerically verify that the dipole mo- 
ment Qi is indeed zero within our numerical error bars. 



3 NUMERICAL INTEGRATION OF THE 
MOND EQUATION 

3.1 Formulation and implementation 

Using spherical spatial coordinates {r, 9,tp}, we consider a 
star represented by a spherically symmetric distribution of 
matter p(r) obtained from a hydrostatic equilibrium model 
in Newtonian theory (polytrope). These models are obtained 
by integrating the hydrostatic spherical equilibrium equa- 
tion (relating the pressure p and the gravitational field) and 
the Newtonian equation for the gravitational field (relat- 
ing the gravitational field and the density distribution p), 
together with a polytropic equation of state (EOS) of the 
form p — Kp~' , where K and 7 are two constants. As it shall 
be shown later in Sec. 3.2, the results do not depend on the 
specific EOS used to obtain this hydrostatic equilibrium, so 
we do not discuss the particular EOS used for obtaining the 
density distribution. It also means that we neglect the ef- 
fect of MOND theory on the structure of the star. This is 
not a severe restriction, since inside the star the gravita- 
tional field is much higher in amplitude than the constant 
ao, making p ~ 1 and Newtonian theory is recovered up to 
very high accuracy. We then solve the MOND equation (1), 
with the boundary conditions given by the constant galactic 
gravitational field g B (see Fig. 1), i.e. lim r ^oo g = g c — g c e. 
In order to be closer to a Poisson-like form of the partial 
differential equation (PDE), we are solving in terms of the 
internal potential u — U — g c ■ x such that lim, — >00 u = 0. 
Defining h = V« = g — g c , the MOND equation (1) becomes 
the following PDE: 

A«=-|- i-^Gp- f^iM gWdthA , (22) 
p(g/a ) { a g J 

with pi = dp/dy being the derivative of the function p with 
respect to its argument y — g/a . 

This PDE is solved numerically, using the library 
lorene which implements spectral methods (for a review in 
the case of numerical relativity see Grandclement & Novak 
2009) in spherical coordinates. In our case (see Fig. 1), the 
fields do not depend on the azimuthal angle <p and the prob- 
lem is axisymmetric. Since Eq. (22) is a non-linear elliptic 
PDE, the algorithm used is the fixed point iteration method 
in which one starts from an initial guess uo(r,9), here the 
solution of the Newtonian Poisson equation Auo = —4iTGp. 
Knowing u n (r,8) at a given iteration step n, one com- 
putes the non-linear source terms in the right-hand-side of 
Eq. (22), say <t(u, p), to obtain a new value of the potential 
u n +i, solving a linear Poisson equation 

Au n+ i = a(u n ,p) , (23) 

and using a relaxation technique 

ttn+l = Xu n +i + (1 - \)u n , (24) 

with A G (0, 1] being the relaxation parameter (often taken 
to be 0.5). This iteration is stopped when the relative varia- 
tion of the potential u between two successive steps becomes 
lower than a given threshold (usually 10~ 12 ). 

The linear equation (23) is solved decomposing each 
field on a truncated base of spherical harmonics y ; m (#, <p) 
(with m — because the problem is axisymmetric) in 
the f?-direction, and Chebyshev polynomials T„(r) in the 
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r-direction. For this last coordinate, a multi-domain tech- 
nique is used, with linear mappings of the r coordinate to 
the interval [— 1, 1], except for the last domain where a map- 
ping of the type s — 1/r allows to treat spatial infinity in 
our computational domain and to impose the right bound- 
ary conditions at r — > oo (see Grandclement & Novak 2009 
for details). Moreover, this multi-domain technique imposes 
that the potential u, as well as its first radial derivative, be 
continuous across the stellar surface. 

Finally, the only modification made to the exact equa- 
tion (22) in order to be numerically integrated is the setting 
fj,(g/ao) = 1 inside the star. Indeed, there are two problems 
here. First, the MOND gravitational field g(r) does not ad- 
mit a regular Taylor expansion near the origin, as does the 
density in Eq. (B6) and the associated Newtonian gravita- 
tional field which follows as 



+o° r 2fc+l 

^^E p^p^ ^tO). (25) 

This can be immediately seen from the fact that in spherical 
symmetry the MOND equation (I) reduces to (j,(g/ao)g = 
<7n- Using (25) we see that the expansion of g when r — > 
starts with a term proportional to the square root of r and 
is therefore not regular; more precisely we have 



9 = 



I + 0(r 2 )] 



(26) 



With our polynomial decomposition of fields in the radial 
direction, where we decompose in Chebyshev polynomials 
T„(r), we therefore get an incompatibility when trying to 
solve the Poisson equation (23) in the vicinity of the origin 
at the centre of the Sun. The second problem is that, at the 
very centre of the star g — > and therefore /i — > 0, which 
makes the division present in the right-hand-side of Eq. (22) 
difficult to handle numerically. 

Here we have chosen to avoid both problems by im- 
posing /j, = 1 in the star so it is entirely Newtonian. This 
approximation may produce a noticeable change in the value 
of n(y) only within a sphere of radius r < 



3an 



47rGpo 



W~ 15 R 



where po is the central density and Rq is the radius of the 
Sun. It is thus supposed to have a completely negligible ef- 
fect on the results. 

The result which is finally used to study the influence 
of the modification of the Newtonian gravity on the orbits 
of planets is the value of the trace-free multipole moments 
Ql defined in Sec. 2. In the case where all the multipole 
moments are induced by the presence of the external field 
g c in the preferred direction e, the situation is axisymmetric 
and all the moments Q l will have their axis pointing in that 
direction e, and we can define the multipole coefficients Qi 
by 



Qi 



Qi e L . 



(27) 



where e denotes the STF part of the product of / unit 
vectors e L — e 11 ■ ■ ■ e H . Then the multipole expansion (14) 



reads as J 

Su(r, 6)=J2 { 2l - 1)!! r ' Ql{r) Pl{c ° s9) ' (28) 

where Pi is the usual Legendre polynomial and is defined 
in Fig. 1. Although from the considerations of Sec. 2 the 
multipole moments Qi should be approximately constant 
within the MOND transition radius ro, here we compute 
them directly from the numerical solution of (1) and shall 
obtain their dependence on r; we shall check numerically 
that Q2(r) and Qa(r) are indeed almost constant in a large 
sphere surrounding the Solar system. With our definition the 
monopole, quadrupole and octupole pieces in the internal 
field are given by 



8u i 



-r Qi cos 6 , 



5u2 = i r 2 Q2 ( cos 2 ' 



8U3 = — i r 3 Q3 cos 6 I cos 3 9 — - cos 



(29a) 
(29b) 

(29c) 



From dimensional arguments, it is possible to see that Qi 
should scale as Qi ~ ao/r 1 ^ 1 [see for instance Eqs. (33) and 
(37)] and therefore 6 



8u ~ ^2 Qi r ' ~ a o r o ^2 ( 
j 1 V 



(30) 



This last expression shows that higher-order multipoles 
should have lower influence on Solar-system planets, for 
which r < ro. This shall be exemplified with the influence 
of the octupole, as compared to that of the quadrupole, in 
Sec. 4.2. 



3.2 Tests of the code 

A first test of our code is the standard comparison between 
the Laplace operator applied to the potential u, and the 
right-hand-side of Eq. (22). In all results displayed here, the 
maximum of this error has always remained lower than 10 -11 
and shall not be considered as an interesting error indicator. 
Next, an indication of the error comes from the use of the 
Gauss theorem: 



g ■ d 2 S = -4ttGM, 



(31) 



where M is the mass of the star, computed from initial data 
at a very high precision. Then, another check (which is how- 
ever not a priori independent) is the asymptotic behaviour 
of the potential u(r, 6) when r — ► 00, which is 



GM 



+ 



(32) 



5 Here, the expansion is denned for any range in radius r, con- 
trary to Sec. 2 and Appendix B, where the multipoles were only 
numbers defined for r — > and not functions of r. 

6 The radial dependence of the contribution of the i-th multipole 
moment is the same as that of the usual MOND spherical effect 
corresponding to a MOND function behaving like /n ~ 1 — ey~ 7 
when y — > 00, with 7 = (I + l)/2 (see the discussion at the end 
of Sec. 1.1). But of course the effect we arc considering here is 
non-spherical. 
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2 

H(y) = y/(l+y ) 



— Test using the Gauss theorem at infinity 
x Test on the asymptotic behavior of u(r,6) 




Figure 2. Relative error of the code on the two tests described in 
Sec. 3.2, for a IMq spherical mass distribution, for the standard 
MOND function fJ.2(y) [Eq. (34b)] and the standard value ao = 
1.2 X 10 _10 m.s- 2 . The abscissa gives the ratio (5) between the 
norm of the asymptotic gravitational field g c and the parameter 
ao. Note that both curves do coincide. 




radius [m] 



Figure 3. Two density profiles for testing the result dependence 
on the equation of state. Both profiles yield a one-solar mass star 
in hydrostatic equilibrium. The two profiles yield very little (same 
level as the overall numerical accuracy ~ 10~ 3 ) difference on our 
results. 



where ^ e = £t(j/ e ) and A c = y e p' e /p e with y c = g c /a (see 
the proof in Appendix A). The two tests are not really in- 
dependent because the Gauss theorem is obtained by inte- 
grating the asymptotic field over a sphere at infinity. The 
main interest of these tests is that while the field is com- 
puted asymptotically on a sphere at infinity, the mass M 
is obtained "locally" by integrating numerically the density 
over the volume of the star. We emphasise that in our code 
the correct asymptotic behaviour (32) comes out directly 
and needs not to be imposed by hands at the beginning of 
the calculation. 

From the use of the mapping s = 1/r (see Sec. 3.1), 
it is numerically possible to check the angular dependence 
of the quantity r x u when r — > oo and thus to estimate 
independently the error on the potential u. Both tests - 
Gauss theorem (31) and asymptotic behaviour (32) — have 
been performed and have given accuracy levels of ~ 10 -3 - 
1CP 4 . Some examples are given in Fig. 2, for a sequence 
of different values for the asymptotic gravitational field g c . 
Both tests give the same numbers, up to five digits. In all 
results shown hereafter, the accuracy level defined by these 
two tests is better than 1CP 2 . 

Among the parameters entering the numerical model, 
we have checked that the details of the density profile p(r) do 
not contribute to the value of the quadrupole moment Q2(r). 
As an illustration, we give here the quadrupole obtained 
with two different density profiles, displayed in Fig. 3, both 
giving a star of exactly one solar mass. With the same other 
parameters [choices of fj,(y), ao and g e ], the relative difference 
in the induced quadrupole is 4 x 1CF 3 , while the error indi- 
cators give a relative numerical uncertainty of ~ 2 x 1CF 3 . 
We therefore consider that the particular form of the density 
profile used to model the star (i.e., the choice of EOS) has 
little influence on the results. 



3.3 Results on the induced quadrupole and 
octupole 

In what follows, unless specified, the mass of the star is that 
of the Sun. As a first result, we show in Fig. 4 the profile 
of the quadrupole induced by the MOND theory through 
the function Qz(r) defined in Eq. (29b). We find that this 
quadrupole is decreasing from the star's neighbourhood to 
zero, on a typical scale of 10000 astronomical units (AU). 
However, this quadrupole can be considered as almost con- 
stant closer to the Sun, as it has a relative variation lower 
than 10~ 4 within 30 AU (see the zoomed region in Fig. 4). 
We shall therefore refer to the quadrupole as a simple num- 
ber, noted Q2(0) or simply Q2, when evaluating its influence 
on the orbits of Solar-system planets. Our numerical results 
for the quadrupole are given in Table 1. 

We now briefly study other multipole moments, namely 
for I — 1 and / = 3. The profiles for the dipole Qi(r) and 
octupole Q3(r), defined by (29a) and (29c), are displayed in 
Fig. 5. We find that near the Sun and the solar planets the 
dipole can be considered as zero, since it is three orders of 
magnitude lower than the typical value for an acceleration in 
the problem (i.e. ao or g e ). Indeed we expect on dimensional 
analysis that Q\ should scale with the MOND acceleration 
ao [see (30)], 

Qi = a qi{rj) , (33) 

where the dimensionless coefficient gi depends on the ratio 
rj between g e and ao as defined by Eq. (5). Our values given 
in Tabic 1 show that gi is actually very small. It means 
that Qi is lower than the numerical error and confirms the 
analytical proof given in Sec. 2 that the dipole moment is 
approximately zero. On the other hand, we find that the 
octupole tends to a non-zero value, because the correspond- 
ing g3, which has also an influence on the orbits of planets 
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Figure 4. Left panel: profile of Q2(r) in the solar system, for a standard choice of function Hi(y) [see Eq. (34a)], ao = 1.2 X 10~ 10 m.s~ 2 
and g e = 1.9 X 10~ 10 m.s~ 2 . The MOND transition radius ro = \fGMjao is shown by a dash-dotted line at r ~ 7100 AU. Right panel: 
zoom of the central region (r ^ 50 AU), where the quadrupole is almost constant (the spikes for r — > arc due to the high value of the 
monopole term, from which it is numerically difficult to extract the quadrupole). Note the difference in the y-axis scales. 




2000 4000 6000 r ° B000 10* 2000 4000 6000 r o 8000 10 4 

Radius r [AU] Radius r [AU] 



Figure 5. Profiles of dipole Q\ (r) (left) and octupole Q3 (r) (right) in the solar system, with same settings as those of Fig. 4. The MOND 
transition radius ro = V GM/ao is shown by a dash-dotted line at r ~ 7100 AU. Numerically we get |Qi(0)| ~9x 10~ 14 m.s~ 2 <C ao 
which shows that the dipole is in fact zero. 



(see the next section), is found to be of the order of one. 
Numerical values of the octupole are also given in Table 1. 

The dependence of the quadrupole moment Q2 upon 
the value of the galactic gravitational field is displayed in 
Fig. 6, for four different coupling functions fi(y). Here we 
consider four cases widely used in the literature: 



general notation for any positive integer n: 



Mi (y) = 



V 



1 + y 



/ioxp(y) = 1 - e~ v , 

/iTcVoS(y) = 



VI + % + 1 



(34a) 
(34b) 
(34c) 
(34d) 



The function /Hi has been shown to yield good fits of galactic 
rotation curves (Famaey & Binney 2005). However because 
of its slow transition to the Newtonian regime it is a prion 
incompatible with Solar System observations. The function 
fj,2 is generally called the "standard" choice and was used in 
fits such as those of Begeman et al. (1991). We here use the 



My) 



y 



(35) 



We include also the function /i exp having an exponentially 
fast transition to the Newtonian regime. The fourth choice 
MToVcS is motivated by the theory TeVeS (Bekenstein 2004). 
One should note that none of these functions, except maybe 
the fourth one, derives from a fundamental physical princi- 
ple. 

One observes that, up to the standard value of g c — 
1.9 x 10~ 10 m.s~ 2 the four curves are quite close, giving an 
interval of values for Q2: 

2.2 x KT 26 s~ 2 < Q 2 < 4.1 x 10~ 26 s~ 2 . (36) 

However, for other choices of the MOND function fi(y), Q2 
can take lower values down to 2.1 x 10~ 27 s~ 2 . In all cases, we 
find that Q2 is positive, which means a prolate deformation 
of the field toward the galactic centre. Two profiles of fi(y) 
seem to give a maximum for Q2 near the standard value 
of g e , namely the standard choice /12 and the exponential 
choice /i exp . For the two other choices of fi(y), we have not 
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Figure 6. Quadrupole in the solar system Q2 = Q'2(0) as a 
function of the external galactic gravitational field g e , for four 
different expressions of the coupling function fJ.(y) [see Eqs. (34) 
for a description]. The standard value of g c = 1.9 X 10~ 10 m.s~ 2 
is shown by a thin dash-double-dotted line. 



Figure 8. Quadrupole moment Q2 as a function of the stellar 
mass M. The M _1 / 2 -law is recovered as expected from Eq. (37a). 



a= 1.2x10 m/s" 



— le-26 - X 



give values of the dipole, quadrupole and octupole near the 
Sun (r 50 AU), where they are observed to be constant, 
for the four different expressions (34) of the interpolating 
function. 

We have also explored different values of ao at fixed 
g c , with the results that Q2 was increasing for small ao as 
a power-law, reaching a maximum value for ao — 1CP 10 — 
10~ 9 m • s~ 2 and then slowly decreasing. Finally, we have 
varied the mass of the central star. Results are displayed 
in Fig. 8 where, in particular, the M~ 1//2 dependence of 
the quadrupole moment is recovered. This is not a surprise 
since on general grounds one expects that the quadrupole 
moment (and similarly the octupole moment) should scale 
like (Milgrom 2009) 7 



O2 = ^92(77): 

n> 



(37a) 
(37b) 



Figure 7. Quadrupole moment Q2 as a function of the index n 
of the MOND function n n {y) defined by Eq. (35). 



been able to increase the value of r\ = g c /ao further than 
10, because the error indicators would become too large and 
the results could not be trusted any longer. Note that in 
Fig. 6, ao has been kept fixed, while we have varied g e . We 
have further explored the dependence of the EFE induced 
quadrupole Q2 on the type of MOND function. 

We have used several different functions of type /i„, as 
defined in Eq. (35), and the results for the variation of Q2 
depending on n are displayed in Fig. 7. From this figure, 
one can notice that the value of Q2 decreases with n, that 
is with a faster transition from the weak-field regime (where 
Mf) ~ V) to the strong field regime {n{y) ~ 1). We have 
been unable to study higher values of n and to determine 
a possible limit for Q2 as n goes to infinity. In Table 1, we 



where ro = y GM/ao is the MOND transition radius and 
02, 03 denote some dimensionless coefficients depending on 
the ratio r\ = g c /ao and on the choice of the interpolating 
function fi. Having obtained the expected dependence on 
the mass in Fig. 8 is another check that our code behaves 
correctly, since we are able to recover the analytic results 
known for this system. The values of the dimensionless mul- 
tipole coefficients 02 and 03, which are expectedly close to 
unity, are reported in Table 1. 



' To compare with the results of Milgrom (2009) for the 
quadrupole, one should note the different conventions in use, 
namely ^ = —Qff and ^Milgrom {rj) = _| g BN (j?) _ Qur rc _ 

suits for the quadrupole are in good agreement with values given 
by Milgrom (2009). 
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Table 1. Numerical values of the dipole Q\, the quadrupole Q2 and the octupole Q3 together with the associated dimensionless quantities 
qi,Q2 and 93 defined by Eqs. (33) and (37). All values are given near the Sun. We use different choices of the function /i(y) defined in 
Eqs. (34) and (35). These values have been obtained for ao = 1.2 X 10~ 1() m • s~ 2 and g c = 1.9 X 10~ 10 m • s~ 2 (and M = 1 M Q ). 



MOND function 


Mife) 






^20 (y) 


Mcxpfe) 


MToVeS (V) 


Ql [ms- 2 ] 


-8.9 x itr 14 


-9.8 x 10~ 14 


-1.1 x 10~ 13 


-1.2 x 10~ 13 


-9.7 x 10~ 14 


-3.5 X 10~ 14 


qi 


-7.4 x 1CT 4 


-8.2 x 10~ 4 


-9.2 x 10~ 4 


-10~ 3 


-8.1 x 10~ 4 


-2.9 x 10~ 4 


Q2 [s- 2 ] 


3.8 x 1(T 26 


2.2 x 10" 26 


7.4 x 10" 27 


2.1 x 10~ 27 


3.0 x 10" 26 


4.1 x 10~ 26 


<?2 


0.33 


0.19 


6.5 x 10~ 2 


1.8 x 10~ 2 


0.26 


0.36 


Qs [m^s- 2 ] 


-1.2 x lO" 40 


-9.3 x 10~ 41 


-4.9 x 10~ 42 


-2.3 x 10~ 42 


-1.2 x lO" 40 


-1.1 x IO- 40 


93 


-1.1 


-0.87 


-4.6 x 10~ 2 


-2.1 x 10~ 2 


-1.1 


-1 



4 EFFECT ON THE DYNAMICS OF SOLAR 
SYSTEM PLANETS 

4.1 Perturbation equations 

In this section, we investigate the consequence for the dy- 
namics of inner planets of the Solar System of the presence 
of an abnormal multipole moment Ql oriented toward the 
direction e of the galactic centre, in the sense of Eq. (27). 
Recall that the domain of validity of this anomaly is ex- 
pected to enclose all the inner Solar System (for distances 
r < ro « 7100 AU), with the quadrupole coefficient being 
constant up to say 50 AU (see Fig. 4). As we have seen, 
the anomaly induces a perturbation on the Newtonian grav- 
itational potential, namely u — un + Su, where the pertur- 
bation function 5u is given for various multipole moments 
by Eqs. (29). Following the standard practice of celestial 
mechanics we denote the perturbation by R = 5u. The 
perturbation function R is such that the perturbing force 
(or, rather, acceleration) acting on the Newtonian motion is 
F = VR. 

The unperturbed Keplerian orbit of a planet around 
the Sun is described by six orbital elements (say ca with 
A = 1, • • • ,6). For these we adopt the semi-major axis a, 
the eccentricity e, the inclination / of the orbital plane, the 
mean anomaly £ defined by £ = n(t — T) where n = 2n/P 
(n is the mean motion, P is the orbital period and T is the 
instant of passage at the perihelion), the argument of the 
perihelion oj (or angular distance from ascending node to 
perihelion), and the longitude of the ascending node Q,. We 
also use the longitude of the perihelion defined by u — lu + Q,. 

The perturbation is a function R(ca) of the orbital ele- 
ments of the unperturbed Keplerian ellipse. With our choice 
for the orbital elements {ca} = {a, e, I,£, cj,f2} the per- 
turbation equations of celestial mechanics read as (see e.g. 
Brouwer & Clemence 1961) 



da 
dt 
de 
dt 
dJ 
dt 
d£ 
dt 

du> 
~dt 



2_dR 

an~d£ ' 

VI - e 2 



a 2 n\/l — e 2 sin / 



; dR_dR 
' d£ du) 

dR dR 
du dQ 
1 - e 2 dR 



cos I- 



2a— + 

da e de 



VI - e 2 dR 
e de 



il 



dR 



i/VT 3 ^ 91 



(38a) 
(38b) 
(38c) 
(38d) 
(38e) 



dQ 
dt 



dR 



72 QI 



(38f) 



a 2 n sin I VT~ 

We shall be interested only in secular effects, so we av- 
erage these equations over one orbital period P; denoting 
the time average by brackets we have 

r 2n 



p J dR 
dt 7T- 

dCA 



1 

2^ 



(1/ 



dR 

dCA 



(39) 



In particular this implies that we shall always have 
(dR/d£) = hence the perturbation equation (38a) gives 
(da/dt) = (at first order in perturbation theory). Indeed, 
a perturbative force of the type F — Vii is conservative, so 
the energy of the orbit is conserved in average and there is 
no secular change in the semi-major axis. 

Let us now apply the perturbation equations (38) to the 
specific case of the perturbation function corresponding to 
the quadrupole anomaly, namely 



R = Su2 = - r 2 



2 « 1 

cos t) — — 



(40) 



x/r being the unit direction 



2 + z 2 , where (x,y,z) are the 



Here cos 9 = e • n, with n = 
of the planet and r 2 = x 2 + 
coordinates of the planet in an absolute Galilean frame cen- 
tred on the focus of the unperturbed Keplerian ellipse, and 
with respect to which the orbital elements {a, e, I, I, w, Q} 
are defined. In the following, to simplify the presentation, 
we shall choose the a;-direction of this Galilean frame to be 
the direction of the galactic centre e — g c /g c . That is, we as- 
sume that the origin of the longitude of the ascending node 
lies in the direction of the galactic centre. This means that 
cos# = x/r so that 

Q 



R=^(2 X 2 -y 2 - Z 2 ) 



(41) 



We express the planet's absolute coordinates (x, y, z) in 
terms of the orbital elements {a, e, I, £, u>, fl} by perform- 
ing as usual three successive frame rotations with angles fi, 
I and to, to arrive at the frame (u, v, w) associated with the 
motion, where (u, v) is in the orbital plane, with u in the 
direction of the perihelion and v oriented in the sense of 
motion at perihelion. The unperturbed coordinates of the 
planet in this frame are 



u — a (cos U — e) , 
v — a\J\ — e 2 sin U , 



(42a) 
(42b) 

w = 0, (42c) 
where U denotes the eccentric anomaly, related to £ by the 
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Kepler equation t = U — esinU. Inserting the resulting ex- 
pression of R into the perturbation equations (38), and per- 
forming the time average (in practice, an average over the 
eccentric anomaly U, after the appropriate change of vari- 
able), then yields 



da 

dt 

de\ 
dt) 



0, 



(43a) 



4n 



cos I cos(2o;) sin(2fi) 



+ sin(2aj) (cos 2 f2 — cos 2 I sin 2 fl)] , (43b) 
Q2 sin I sin Q 



_ [(2 + 3e 2 + 5e 2 cos(2w)) cos Q 



4n\/T— e 

— 5e 2 cos/sin(2u;) sinfi] , 



(43c) 



-) = n+ ^"{-15(1 + e 2 )cos(2a;) [2-2 cos(27) 



d£\ 

dt) ' 96n 

+ cos(2J - 2Q) + 6cos(2n) + cos(2/ + 2Q)\ 
+ 12 (-8 + 3e 2 ) cos 2 ft 
-12 (1 - 6e 2 + (7 + 3e 2 ) cos(27)) sin 2 ft 
+20 (2 - 3e 2 

+6(1 + e 2 ) cos I sin(2a;) sin(2ft)) } , (43d) 



/da; 
\ dt 



[2 - 2e 2 



dft 
dt 



4n\A — e 2 

+ (-l + e 2 ) (3 + 5cos(2u;))cos 2 ft 
5 2 

+ — e cos J sin(2cj) sin(2ft) 



+5(1 - e 2 )cos/sin(2w)sin(2ft)] , (43e) 
Q 2 sin ft 



[5e 2 cosftsin(2w) 



4n%/l — e 2 

+ cos J (-2 - 3e 2 + 5e 2 cos(2cj)) sin ft] (43f) 



These equations are general and give in particular the pre- 
cession of the node (dft/dt) and the perihelion precession 
(dcj/dt) or, rather, (dCo/dt) where Co — u> + ft is the longi- 
tude of the perihelion. 

In order to make some estimate of the magnitude of 
the effect, let us approximate the direction of the galactic 
centre (which is only 5.5 degrees off the plane of the ecliptic) 8 
as being located in the plane of the orbit; consequently we 
choose I = 0. In this case Co = 10 + ft is the relevant angle 
for the argument of the perihelion. The non-zero evolution 
equations then become 



Ide 
\dt 



d£ 
dt 

do; 
717 



An 

Q 



sin(2a3) , 



(44a) 



= n - f^- [7 + 3e 2 + 15(1 + e 2 ) cos(2cD)] , (44b) 



Q2VT 



4n 



1 + 5 cos 



(2u>)] . 



(44c) 



We recall that with our notation uj is the azimuthal an- 
gle between the direction of the perihelion and that of the 



galactic centre (approximated to lie in the orbital plane) . Of 
particular interest is the secular precession of the perihelion 
(do;/di) due to the quadrupole effect which we henceforth 
denote by 9 



Q2VT 



4n 



1 + 5cos(2oi)j 



(45) 



The precession is non-spherical, in the sense that it depends 
on the orientation of the orbit relative to the galactic centre 
through its dependence upon the perihelion's longitude uj. 
The effect scales with the inverse of the orbital frequency 
n = 2-k/P and therefore becomes more important for outer 
planets like Saturn than for inner planets like Mercury. This 
is in agreement with the fact that the quadrupole effect we 
are considering increases with the distance to the Sun (but 
of course will fall down when r becomes appreciably compa- 
rable to ro, see Fig. 4). 

We have also computed the secular planetary precession 
induced by the octupole moment Q3, for which the pertur- 
bation function reads with the same conventions 

R = -i r 3 Q 3 (cos 3 6 - I cosfl) = -^x (2x 2 -3y 2 -3z 2 ) . 

(46) 

Redoing the perturbation analysis we find the octupolar pre- 
cession in the case I = as 



A3 



32en 



[2 - 13e 2 + 35e 2 cos(2cD)] cos Co . (47) 



4.2 Numerical evaluation of the planetary 
precession 

We now compute numerically this effect for the various plan- 
ets of the Solar System; the relevant orbital elements for the 
planets we used in this calculation are provided in Table 2. 

Our numerical values for the quadrupole and octupole 
anomalies A2 and A3 are reported in Tables 3 and 4 respec- 
tively. As we see from Table 3 the quadrupolar precession 
A2 is in the range of the milli-arc-second per century which 
is not negligible. In particular it becomes interestingly large 
for the outer gaseous planets of the Solar System, essen- 
tially Saturn, Uranus and Neptune. The dependence on the 
choice of the MOND function /i is noticeable only for func- 
tions (i n (y) defined by (35) with large values of n, where the 
effect decreases by a factor ~ 10 between n — 2 and n = 20. 
However, the octupolar precession A3 in Table 4 is much 
smaller, being rather in the range of the micro-arc-second 
per century. 

We give for comparison in Table 5 the best published 
postfit residuals for any possible supplementary precession 
of planetary orbits (after the relativistic precession has been 
duly taken into account), which have been obtained from 
global fits of the Solar System dynamics by Pitjeva (2005) 
and Fienga et al. (2009). In particular the INPOP planetary 
ephemerides by Fienga et al. (2009) use information from the 
combination of very accurate tracking data of spacecrafts 
orbiting different planets. 

From Table 5 we see that our results for A2 are smaller 



8 The latitude /3 and longitude A of the galactic centre in the stan- 
dard ecliptic coordinate system are f3 = —5.5° and A = —93.2° 
(see e.g. Allen 1999). 



9 The result found by Milgrom (2009) for this effect [his 
Eqs. (37)-(38)] looks more complicated, but can be simplified 
and is seen to be equivalent to ours. 
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Table 2. Orbital elements of planets used in the computation of the abnormal precession rates A 2 and A3. The longitude of the perihelion 
u> is defined here with respect to the direction of the external galactic field (assuming that the galactic centre lies within the plane of the 
ecliptic). Thus we pose ui = u}\n en — A, where A = —93.2° is the longitude of the galactic centre and a) Alien is given in Allen (1999). 







Mercury 


Venus 


Earth 


Mars 


Jupiter 


Saturn 


Uranus 


Neptune 


a 


[AU] 


0.397 


0.723 


1.00 


1.523 


5.203 


9.537 


19.229 


30.069 


e 




0.206 


0.007 


0.017 


0.093 


0.048 


0.054 


0.047 


0.009 


P 


M 


0.24 


0.62 


1.00 


1.88 


11.86 


29.46 


84.01 


164.8 




[deg] 


170.7 


224.7 


196.1 


69.2 


108.0 


185.6 


264.2 


138.2 



Table 3. Results for the precession rates of planets A2 due to the quadrupolc coefficient Q2- We use the values for Q2 corresponding to 
various MOND functions defined in Eqs. (34)-(35) (see Table 1). All results are given in miZB-arc-seconds per century (mas/cy). 



Quadrupolar precession rate A2 in mas/cy 



MOND function 


Mercury 


Venus 


Earth 


Mars 


Jupiter 


Saturn 


Uranus 


Neptune 




0.04 


0.02 


0.16 


-0.16 


-1.12 


5.39 


-10.14 


7.93 


M2(y) 


0.02 


0.01 


0.09 


-0.09 


-0.65 


3.12 


-5.87 


4.59 


Ms (?/) 


7 x 10~ 3 


3 x 10~ 3 


0.03 


-0.03 


-0.22 


1.05 


-1.97 


1.54 


M20 (y) 


2 x 10~ 3 


10~ 3 


9 x 10~ 3 


-9 x 10~ 3 


-0.06 


0.3 


-0.56 


0.44 




0.03 


0.02 


0.13 


-0.13 


-0.88 


4.25 


-8.01 


6.26 


MTeVeS (y) 


0.05 


0.02 


0.17 


-0.17 


-1.21 


5.81 


-10.94 


8.56 



or much smaller than the published residuals except for 
the planets Earth, Mars and Saturn. The case of Saturn 
is interesting because the INPOP planetary ephemerides by 
Fienga et al. (2009) publish an offset with respect to the gen- 
eral relativistic value for the precession: namely they quote 
— 10 ± 8 mas/cy as obtained from the Cassini tracking data 
and 200 ± 160 mas/cy from the VEX data. The values we 
find for A2 are smaller but grossly within the range of these 
postfit residuals. However we find that A2 is positive for 
Saturn while the offset reported in Fienga et al. (2009) from 
Cassini data is negative. But note that this may not be rel- 
evant because the INPOP ephemerides are used by Fienga 
et al. (2009) to detect the presence of an eventual abnormal 
precession, not to adjust the value of that precession. 

In the case of the Earth the INPOP ephemerides find a 
tight constraint ± 0.016 coming from the Jupiter VLBI 
data (Fienga et al. 2009). This constraint seems already 
to exclude most of our obtained values for A2, except for 
MOND functions of type fj, n (35) with rather large values 
of n. In particular, one needs to take n > 8 in order to 
have an EFE precession compatible with this constraint on 
the Earth orbit. On the other hand we note that for all the 
planets the octupolar precession rate A3 is very small and 
clearly completely negligible in current fits of the Solar Sys- 
tem dynamics. 

Thus it seems that in the case of the Earth the con- 
straint from the Jupiter VLBI data is already quite tight, 
and it is difficult to accommodate our anomalous quadrupo- 
lar precession rate A2 for most choices of MOND function 
n(y). However let us remark that the postfit residuals in 
Table 5 are obtained by adding by hands an excess of pre- 
cession for the planets and looking for the tolerance of the 
data on this excess. But in order to really test the anoma- 
lous quadrupolar precession rate A2, one should consistently 
work in a MOND picture, i.e. consider also the other effects 
predicted by this theory, like the precession of the nodes, 



the variation of the eccentricity and the inclination, and so 
on. Then one should perform a global fit of all these effects 
to the data; it is likely that in this way the quantitative 
conclusions would be different. 

On the other hand, as we have commented in the Intro- 
duction, it is non-trivial to know if testing MOND-like the- 
ories (including TeVeS and the model proposed by Blanchet 
& Le Tiec 2008, 2009) in the Solar System could invalidate 
those theories, if they represent only some phenomenologi- 
cal models describing the weak field regime of gravity and 
as such should not be extrapolated as they are in the strong 
field of the Solar System. 



5 CONCLUSION 

In MONDian gravity, the influence of the external galac- 
tic gravitational field onto the orbits of Solar-system plan- 
ets appears mainly through the presence of a quadrupolar 
perturbation to the gravitational potential, as it has been 
shown through a multipolar analysis of the solution of the 
MOND equation. This contribution has been computed us- 
ing a spectral, highly accurate, numerical Poisson solver to 
obtain a solution of the equation (1). This numerical solu- 
tion has been extensively tested and compared to known an- 
alytical formulas, and all the expected properties have been 
recovered. Using the most commonly used MOND coupling 
functions jj,(y) and value of the MOND acceleration ao, wc 
have obtained quadrupole contributions that are compati- 
ble with those published by Milgrom (2009). We have also 
extended this calculation to include octupole contributions. 
With the derivation of detailed perturbation equations, the 
influence of the quadrupole and octupole on the perihelion 
precession of planets have been quantified and compared to 
two sets of observations, in particular the INPOP planetary 
ephemerides by Fienga et al. (2009). 
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Table 4. Results for the precession rates of planets A3 due to the octupole coefficient Q3. Beware of the results given here in micro- 
arc-seconds per century (/ias/cy). 



Octupolar precession rate A3 in /ias/cy 



MOND function 


Mercury 


Venus 


Earth 


Mars 


Jupiter 


Saturn 


Uranus 


Neptune 




2 X 1(T 3 


0.2 


0.2 


-0.03 


1.4 


19.5 


12.1 


1499 


M2(y) 


2 x itr 3 


0.1 


0.2 


-0.03 


1.1 


15.1 


9.4 


1162 


(y) 


10- 4 


5 x 10~ 3 


0.01 


-2 x 10~ 3 


0.06 


0.8 


0.5 


61.2 




5 x icr 5 


2 x 10~ 3 


5 x 10~ 3 


-7 x 10~ 4 


0.03 


0.4 


0.2 


28.7 


Moxp (y) 


2 x itr 3 


0.2 


0.2 


-0.03 


1.4 


19.5 


12.1 


1499 


MToVeS (.y) 


2 x 10~ 3 


0.2 


0.2 


-0.03 


1.3 


17.9 


11.1 


1374 



Table 5. We reproduce some published residuals for any supplementary orbital planetary precession. All results are given in milli-arc- 
seconds per century (mas/cy). 







Postfit residuals for A = (dui/dt) 


in mas/cy 








Origin 


Mercury 


Venus Earth Mars 


Jupiter 


Saturn 


Uranus 


Neptune 


Pitjeva (2005) 
Ficnga et al. (2009) 


-3.6 ±5 
-10 ±30 


-0.4 ±0.5 -0.2 ±0.4 0.1 ±0.5 
-4 ±6 0± 0.016 0±0.2 


142 ± 156 


-6±2 
-10 ±8 


± 2 X 10 4 


± 2 X 10 4 



At first sight, these observational results seem to ex- 
clude the presence of a MOND-EFE induced quadrupolc, 
for some choices of the MOND function /x and scale do- This 
is particularly true for the Earth orbit, where constraints 
from the Jupiter VLBI observations are very strong. How- 
ever, for other choices of the MOND function the observa- 
tional constraints on the precession of the Earth orbit are 
still compatible with the computed effect. These compatible 
MOND functions share the property of having a very sharp 
transition from the weak-field (modified) to the strong field 
(Newtonian) regime. Moreover one should be very cautious 
and keep in mind that: 



1. The observational constraints on the Solar-system or- 
bits have been obtained by a global fit, using a fully relativis- 
tic first-post-Newtonian model. Taking only into account a 
particular MOND effect like precession without considering 
a fully-MONDian model is in principle incoherent since all 
the other effects should also be considered; 

2. The MOND theory is a phenomenological approach try- 
ing to describe gravity in the weak field regime (where evi- 
dence for dark matter arises) . It is therefore unclear if such 
theory can be applied to stronger gravitational fields and, in 
particular, in the Solar System. Conversely, the constraints 
obtained here in strong field may not be relevant for the 
weak-field case unless some more fundamental theory, valid 
in both regimes, is known. 



In any case, further studies are to be done if one wants 
to obtain more stringent conclusions about constraints im- 
posed by Solar-system observations onto MOND-like theo- 
ries. More precise observations could also give valuable in- 
formations about an eventual EFE due to MOND theory 
and also restrict the number of possible MOND functions 
that are compatible with the observations. 
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APPENDIX A: ASYMPTOTIC BEHAVIOUR OF 
THE FIELD 

In this Appendix we prove that the asymptotic behaviour of 
the solution of the Bekenstein-Milgrom equation (1) at large 
distances from an isolated matter source in the presence of 
an external field g a , is given by Eq. (32). The proof has 
already been given by Bekenstein & Milgrom (1984); here 
we present a slightly different derivation. For simplicity we 
suppose that the source is spherically symmetric, so that the 
asymptotic form of the MOND potential is axisymmetric 
around the direction of the external field and of the type 



U 



r V f* 



(Al) 



where f(6) is a function of the azimuthal angle 9 measured 
from the direction e of the external field. We naturally as- 
sume that the function f(9) is regular. With the ansatz (Al) 
we obtain successively the leading order terms in the expan- 
sion of various quantities as 

sin 6 



cos9 

-J + zrrnJ» 



+ o 



p, = H e 1 + 



MP = Mc 9. 



g e r-> 



+ 1 



cos<9 



cos Of — sin 9fe 



—X e f cos 9 — X e fe sin 9 - 



]) + °(*)- 



(A2a) 



fe 
sinO 



(A2b) 



(A2c) 



where we denote fe = df/d9, p e = p(y c ), A c = y c p' e /p e 
(with j/ c = g c /ao and p' c — dp(y c )/dy c ). Next, we impose 
that di(fig l ) = holds asymptotically far from a localised 
matter distribution. This, together with the fact that / is 
a regular function of the azimuthal angle 9, yields the first- 
order differential equation 

(1 + Ac sin 2 9) fe + A c sin 9 cos 9f = , (A3) 

whose solution is given in terms of a constant K by 

/ = ; = 5 • (A4) 



VT+ A c sin 2 9 

Finally we determine this constant by using the Gauss the- 
orem (31), i.e. 



j pg i r 2 n i dO, = -4ttGM, 



(A5) 



where we choose a coordinate sphere since the result does 
not depend on the chosen surface at infinity. Inserting into 
(A5) the asymptotic expansion of the field given by (A2c), 
together with the solution (A4) found for /, we obtain 

K=™, (A6) 

He 



which establishes the looked-for formula (32). 



APPENDIX B: EXPRESSION OF THE 
MULTIPOLE EXPANSION 

The result (19) can also be proved directly, without resorting 
to the computation of surface integrals of Sec. 2. We first 
expand 8u — u — iin when r —¥ using Taylor's formula as 



(Bl) 



Here the derivative operator 8l is a priori non STF but we 
can decompose it into a sum of STF components using the 
formula (see the Appendix A of Blanchet & Damour 1986 
for a compendium of useful formulas) 



[1/2] 



9l = a{ Spa d L -2K) A fc . 



(B2) 



where k ranges from to the integer part of 1/2 denoted [1/2], 
where b~2K denotes a product of k Kronecker symbols carry- 
ing 2k indices chosen among those of L = i\ ■ ■ ■ ij, where the 
parenthesis indicate full symmetrisation over the I indices 
i\ ■ ■ • ii , where A = da denotes the ordinary Laplacian, and 
where the coefficients depending on I and k are given by 

, l\ (2Z-4fc + l)H 

k {l-2k)\ (2Jfc)!!(2/-2fc + l)!! ' { ' 

Using now the fact that ASu — in a neighbourhood of the 
origin (r < e) we see that only the term k — will survive 
in the latter sum, and since a' k=0 = 1 we conclude that we 
can replace in fact the derivative operator in (<9l<5u)(0) by 
its STF projection Bl, hence 



8u = Y,J\ xL ^ l5u )(°)- 



(B4) 



Next, we prove that the contribution from the Newtonian 
potential «n in (B4) is purely monopolar. This comes from 
our assumption that the density of ordinary matter p in the 
Sun is both regular (i.e. C°°) and spherically symmetric, 
i.e. depending only on r (neglecting the back-reaction of the 
field on its source). As we shall now prove, this means that 
the Taylor expansion of p(r) in spherical symmetry must 
contain only even powers of r. The density being regular its 
expansion when r — ¥ is given by Taylor's formula as 



(B5) 



Next we again use the formula (B2) to decompose (B5) into 
a sum of STF pieces. Because p is spherically symmetric 
all these STF pieces must be zero except the monopolar 
contribution; hence only the terms for which / = 2k survive 
in the latter decomposition. Using the value of the coefficient 
l/(2k + 1) we are then left with 



2k 
a k 



E 



(2ft + 1) 



(A k p)(0). 



(B6) 



As we see p(r) admits an expansion containing only even 
powers of r. Now the spherically-symmetric Newtonian po- 
tential «n will also admit an expansion only in even powers 
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of r, so we write accordingly wn(t') = f(r 2 ). Acting on a 
function of r 2 the STF derivative operator gives dhfif 2 ) = 
n L (2r)' f( l \r 2 ) where denotes the l-th derivative of this 
function (see Blanchet & Damour 1986). Hence we find that 
Olun is a perfectly regular function and we conclude that 
(9l«n)(0) = except for the monopole case I = 0. We can 
therefore replace (§lSu)(0) in (B4) by (8lu)(0) for any I > 1 
and we have finally established the result 

+ °o i 

Su = -u N (0) + Y,j7X L (dLu){0). (B7) 
i=o L 



which is the same as in Eq. (19). 



